/*
 * @file random.h
 *
 * Accumulative header for random number generation
 *
 * @author Dahua Lin
 */

#ifndef _MSC_VER
#pragma once
#endif

#ifndef GMLPP_RANDOM_H_
#define GMLPP_RANDOM_H_

#include <gmlpp/base/arma.h>

#include <gmlpp/core/sfmt.h>
#include <gmlpp/core/basic_rng.h>
#include <gmlpp/core/discrete_sampler.h>

namespace gml
{
	// default random number generator

	typedef sfmt_19937 default_rstream_t;
	typedef RandomGen<default_rstream_t> default_rng_t;


	/******************************************************
	 *
	 *  Random number generation with matrix/vector
	 *
	 ******************************************************/

	// real uniform distribution

	template<class Rng>
	inline void fill_randu(Rng& rg, mat& X)  // U(0, 1)
	{
		const uword n = X.n_elem;
		double *x = X.memptr();
		for (uword i = 0; i < n; ++i)
			x[i] = rg.randu();
	}

	template<class Rng>
	inline void fill_randu(Rng& rg, mat& X, const double b) // U(0, b)
	{
		const uword n = X.n_elem;
		double *x = X.memptr();
		for (uword i = 0; i < n; ++i)
			x[i] = rg.randu() * b;
	}

	template<class Rng>
	inline void fill_randu(Rng& rg, mat& X, const double a, const double b) // U(a, b)
	{
		const uword n = X.n_elem;
		double *x = X.memptr();

		const double d = b - a;
		for (uword i = 0; i < n; ++i)
			x[i] = rg.randu() * d + a;
	}

	// normal distribution

	template<class Rng>
	inline void fill_randn(Rng& rg, mat& X) // N(0, 1)
	{
		const uword n = X.n_elem;
		double *x = X.memptr();
		for (uword i = 0; i < n; ++i)
			x[i] = rg.randn();
	}

	template<class Rng>
	inline void fill_randn(Rng& rg, mat& X, const double sigma) // N(0, sigma^2)
	{
		const uword n = X.n_elem;
		double *x = X.memptr();
		for (uword i = 0; i < n; ++i)
			x[i] = rg.randn() * sigma;
	}

	template<class Rng>
	inline void fill_randn(Rng& rg, mat& X, const double mu, const double sigma) // N(mu, sigma^2)
	{
		const uword n = X.n_elem;
		double *x = X.memptr();
		for (uword i = 0; i < n; ++i)
			x[i] = rg.randn() * sigma + mu;
	}

	// exponential distribution

	template<class Rng>
	inline void fill_rande(Rng& rg, mat& X) // Exp(1)
	{
		const uword n = X.n_elem;
		double *x = X.memptr();
		for (uword i = 0; i < n; ++i)
			x[i] = rg.rande();
	}

	template<class Rng>
	inline void fill_rande(Rng& rg, mat& X, const double scale) // Exp(1/scale)
	{
		const uword n = X.n_elem;
		double *x = X.memptr();
		for (uword i = 0; i < n; ++i)
			x[i] = rg.rande() * scale;
	}

	// Gamma distribution

	template<class Rng>
	inline void fill_randg(Rng& rg, mat& X, const double alpha) // Gamma(alpha, 1)
	{
		const uword n = X.n_elem;
		double *x = X.memptr();
		for (uword i = 0; i < n; ++i)
			x[i] = rg.randg(alpha);
	}

	template<class Rng>
	inline void fill_randg(Rng& rg, mat& X, const double alpha, const double scale) // Gamma(alpha, 1/scale)
	{
		const uword n = X.n_elem;
		double *x = X.memptr();
		for (uword i = 0; i < n; ++i)
			x[i] = rg.randg(alpha) * scale;
	}


	/******************************************************
	 *
	 *  Random number generation with matrix/vector
	 *
	 ******************************************************/

	template<class Rng, typename T>
	void rand_suffle(Rng& rg, Mat<T>& x)
	{
		const uword n = x.n_elem;
		const uword n_minus_1 = n - 1;

		T *xp = x.memptr();

		for (uword i = 0; i < n_minus_1; ++i)
		{
			uword j = i + rg.randi(n - i);

			T tv = xp[i];
			xp[i] = xp[j];
			xp[j] = tv;
		}
	}

	template<class Rng>
	inline Row<uword> rand_suffle_row(Rng& rg, uword n)
	{
		Row<uword> x(n);
		for (uword i = 0; i < n; ++i) x[i] = i;
		rand_suffle(rg, x);
		return x;
	}

	template<class Rng>
	inline Col<uword> rand_suffle_col(Rng& rg, uword n)
	{
		Col<uword> x(n);
		for (uword i = 0; i < n; ++i) x[i] = i;
		rand_suffle(rg, x);
		return x;
	}


}

#endif /* RANDOM_H_ */
